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Abstract. A particle scheme for scalar conservation laws in one space 
dimension is presented. Particles representing the solution are moved 
according to their characteristic velocities. Particle interaction is re- 
solved locally, satisfying exact conservation of area. Shocks stay sharp 
and propagate at correct speeds, while rarefaction waves are created 
where appropriate. The method is variation diminishing, entropy de- 
creasing, exactly conservative, and has no numerical dissipation away 
from shocks. Solutions, including the location of shocks, are approxi- 
mated with second order accuracy. Source terms can be included. The 
method is compared to CLAWPACK in various examples, and found to 
yield a comparable or better accuracy for similar resolutions. 



1. Introduction 

Conservation laws are important models for the evolution of continuum 
quantities, describing shocks and rarefaction behavior. Fundamental math- 
ematical properties are global and local conservation, the presence of sim- 
ilarity solutions, and the method of characteristics. Successful numerical 
methods employ these properties to their advantage: Finite difference meth- 
ods yield correct shock speeds if applied in conservation form. Finite vol- 
ume methods are fundamentally based on conservation properties. Godunov 
schemes [10], front tracking methods [12], and many related approaches, ap- 
proximate the global solution by local similarity solutions. The method of 
characteristics is used in the CIR method [4] in combination with an inter- 
polation scheme. Although for scalar equations it provides a direct formula 
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for the solution (where it is smooth), it is less popular, since it does not 
possess conservation properties. Consequently, basic CIR schemes do not 
yield correct shock speeds. 

Many commonly used numerical methods operate on a fixed Eulerian 
grid. Advantages are simple data structures and an easy generalization to 
higher space dimensions. Eulerian schemes can be constructed by tracking 
the "correct" approximate solution for a short time step, either by solving 
local Riemann problems (Godunov [10]) or by tracing characteristics (CIR), 
followed by an interpolation step, at which the solution is remapped onto the 
fixed grid. This "remeshing" step generally yields numerical dissipation and 
dispersion. Since the shortest interaction time between shocks or character- 
istics determines the global time step, remeshing is performed unnecessarily 
in many places. In practice, Eulerian methods require sophisticated schemes 
to obtain solutions with sharp features, but without creating oscillation. Fi- 
nite volume methods are equipped with limiters ^20j, while finite difference 
methods use nonlinear approximations, such as ENO [11] or WENO [16] . 

An alternative approach is to abandon the Eulerian property, and thus 
avoid remeshing. Godunov methods become front tracking methods, at 
least in one space dimension. While in the former the interaction of shocks 
is avoided by remeshing, in the latter it is resolved after approximating 
the flux function by a piecewise linear function. By construction, front 
tracking is successful when representing shocks, but cumbersome when ap- 
proximating smooth parts of the solution. Similarly, CIR methods become 
Lagrangian particle methods. Particles carry function values and move with 
their characteristic velocities. As motivated in [8j, this provides a simple and 
accurate solution method for conservation laws, without ever approximat- 
ing derivatives. However, particle management is required, for two reasons: 
First, neighboring particles may depart from each other, resulting in poorly 
resolved regions. This is prevented by inserting particles into gaps. Sec- 
ond, particles may collide. If left unchecked, such a shock event leads to 
a "breaking wave" solution. This is prevented by merging particles upon 
collision. 

Lagrangian particle methods have been successfully applied in the sim- 
ulation of fluid flows. Examples are vortex methods [2], smoothed particle 
hydrodynamics (SPH) [17^ [9} [18]. or generalized SPH methods [6]. The so- 
lution is approximated on a cloud of points which move with the flow, thus 
the governing equations are discretized in their more natural Lagrangian 
frame of reference. In specific applications, more accurate solutions may be 
obtained than with fixed grid approaches. In addition, with particles local 
adaptivity is a straightforward extension. 

The particle method presented here combines the method of character- 
istics (where the solution is smooth) and particle merges (at shocks). The 
evolution of area between neighboring particles is derived from local simi- 
larity solutions. The method is designed to conserve area exactly. 
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1.1. Formulation of the Particle Method. The simplest form of a one 
dimensional scalar conservation law is 

ut + f{u)x = 0, u{x,0)=uo{x) (1) 

with /' continuous. The characteristic equations [7] 

\u = ^ ' 

yield the solution (while it is smooth) forward in time: At each point 
(xo, tio(xo)) a characteristic curve x{t) = xo + f'{uo{xo))t starts, carrying the 
function value u{x{t),t) = uo{xo). While the particle method is presented 
here for the simple case ([1]), the method of characteristics applies in more 
general cases, such as space-dependent flux functions and source terms (see 
Sect. [7]). When characteristic curves collide, a shock arises. It moves at a 
speed so that area (under the function u{-,t)) evolves correctly with respect 
to ([1]). The Rankine-Hugoniot condition [7] follows from this principle. If 
the flux function / is convex or concave between the left and right state 
of a discontinuity, then the solution forms either a shock or a rarefaction 
wave, i.e. a continuous function connecting the two states. Otherwise, com- 
binations of shocks and rarefactions can result. These physical solutions are 
defined by a weak formulation of ([T]) accompanied by an entropy condition 

The first step in a particle method is to approximate the initial function uq 
by a finite number of points xi < • • • < Xm with function values ui, . . . , Um- 
In Sect.m we present strategies on how to sample the initial function "well". 
The evolution of the solution is found by moving each particle Xi with speed 
f'{ui). This is possible as long as there are no "collisions" between particles. 
Two neighboring particles Xi{t) and Xj+i(t) collide at time t + Ati, where 

f'[Ui+l) - f'{Ui) 

A positive Ati indicates that the two particles at Xi and Xj+i will eventually 
collide. Thus, t + Atg is the time of the next particle collision, where 

Ats = min {{Ati\Ati > 0} U oo} . (4) 

For any time increment At < At^ the particles can be moved directly to 
their new positions Xi + f'{ui)At. Thus, we can step forward in time an 
amount At^- Then, at least one particle will share its position with another. 
To proceed further, we merge each such pair of particles. If the collision 
time Ati is negative, the particles depart from each other. Although at 
each of the particles the correct function value is preserved, after some time 
their distance may be unsatisfyingly large, as the amount of error introduced 
during a merge grows with the size of the neighboring gaps. To avoid this, we 
insert new particles into large gaps (see Sect. 13. ip before merging particles. 

In this paper, we present a method of merging and inserting particles in 
such a way that shocks move at correct speeds, and rarefactions have the 
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correct shape. The strategy is based on mnnicking the evolution of area for 
a conservation law, as is derived in Sect. [2j The definition of an area func- 
tion gives rise to a natural interpolation between neighboring Lagrangian 
particles. As presented in Sect. [3l particle management can then be done to 
conserve area exactly. The resulting particle method is shown to be TVD. 
Since the characteristic equation is solved exactly, and particle management 
is purely local, the method yields no numerical dissipation (where solutions 
are smooth) and correct shock speeds (where they are not). Specific strate- 
gies for sampling the initial data are discussed in Sect. HI 

In the remaining sections, the method is analyzed and generalized. In 
Sect. [5l we prove that the numerical solutions satisfy the Kruzkov entropy 
condition, thus showing that the method yields entropy solutions for convex 
entropy functions. In Sect. [6] we present how non-convex flux functions can 
be treated. Strategies to include sources are presented in Sect. [71 In Sect. [51 
we apply the method to examples and compare it to traditional finite volume 
methods using CLAWPACK [3]. Conclusions are drawn in Sect. [H as well 
as possible applications and extensions of the method outlined. 

2. Evolution of Area for Scalar Conservation Laws 
Consider a one dimensional scalar conservation law 

ut + f{x,u)x = 0, u{x,0) = uo{x) . (5) 
Its characteristic equations [7] 

\x = fuix,u) 
[U = -fx{x,U) 

yield the movement and change of function value of a particle. Let u{x,t) 
be a solution of ([5]). The change of area between two fixed points xi and X2 
is solely given by the flux function / as 



d 



X2 



— j u{x,t)dx = f{xi,u{xi,t)) - f{x2,u{x2,t)) = -[f]ll . (7) 



uyu., L)KX^ — J 1^X1, LLyj^i, L)) — J yu.2,uyd.2, i- J J — — IJ 

In contrast, the change of area between two Lagrangian particles (xi(t), ui{t)) 
and {x2{t),U2{t)), i.e. points that move according to ([6]), is given by 

d /■^2{t) 

— - / U{x,t)dx = {fu{x2,U2)u2 - f{x2,U2)) - {fu{xi,Ui)ui - f{xi,Ul)) 

= F{x2,U2)-Fix,,u,) = [F][2%f^ , (8) 

where F = fuU — f is the Legendre transform of /. That is, / is a Hamilton- 
ian of the dynamics ([6]), and -F is a Lagrangian. Equation ^ (respectively 
dSl)) yields the change of area between two Eulerian (Lagrangian) points, 
only by knowing the flux / (the Lagrangian F) at the two points. Hence, 
in the same fashion as ([7]) can be used to construct a conservative fixed grid 
method, we use ([8]) to construct a conservative particle method. 
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Consider an area value Ai(t) associated with each particle, such that 
[j4]i^^+^ = Ai+i — Ai is the area between Xi and Xj+i. Assume the values Ai 
are known at t = 0. Then we can find the areas at any time by solving the 
system arising from equations ([6]) and ([8]) 

Ui = -fx{xi,Ui) (9) 
Ai = F{xi,Ui) . 

Remark 1. While / = (since / is a Hamiltonian of the dynamics), in 
general F ^ 0. However, if the flux function satisfies 

fxufuU fuufxU fxfu — ) (-^0) 

then F = (by the chain rule). Property (llOp is satisfied for instance if 
/ = f{u) or /(x, u) = ip{x)u^ . If = 0, the evolution of area is particularly 
simple, namely Ai changes at a constant rate -Fj. 

2.1. Space-independent Flux. Henceforth we only consider flux func- 
tions that are independent of the spatial variable, / = f{u). Thus, by 
Rem. [H the area between two Lagrangian points changes linearly, as does 
the distance between them 

- u{x,t)dx = [Fiu)]ll , (11) 

"I Jxi{t) 

f^{x2{t) - X,{t)) = X2{t) - Xl{t) = f'M - fim) = [nu)]2 . (12) 

If the two points xi and X2 move at different speeds, then there is a time 
to (which may be larger or smaller than t) at which they have the same 
position. This assumes that they remain characteristic points between t 
and to, i.e. they do not interact with shocks. At time to, the distance and 
the area between the two points vanish. From (jlip and (|12p we have that 

r'^'\{x,t)dx = {t-to)-[F{u)]ll , 

Jxi{t) 

X2{t)-x,{t) = {t-to)-[f'{u)]ll . 
In short, the area between two Lagrangian points can be written as 

/ u{x,t)dx = {x2{t) - xi{t))af{ui,U2) , (13) 

Jxi{t) 

where af{ui,U2) is the nonlinear average function 

If there is only one flux function, we drop the subscript, and simply write 
a(Mi, U2). The integral form shows that a is indeed an average of u, weighted 
by /". The evolution of area (fTSll is independent of the specific solution. 
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since by assumption we have excluded all solutions for which a shock would 
interact with either characteristic point. The following lemma describes 
some properties of the nonlinear average a(-, ■)• 

Lemma 2. Let f be strictly convex in [Uj^,u^], that is, f" > in {u^,u^). 
Then for all ui,U2 G [u^,u^], the average is. . . 

(1) the same for f and —f; 

(2) symmetric, a{ui,U2) = a{u2,ui); 

(3) an average, i.e. a{ui,U2) G {ui,U2), for ui ^ U2; 

(4) strictly increasing in both ui and U2 ; and 

(5) continuous at ui = U2, with a{u,u) = u. 

Due to the first two properties, we can assume WLOG that /" > and 
ui < U2 whenever convenient. 

Proof. We prove the claims in turn. 
(1,2) Multiplying both numerator and denominator by —1 yields the proof: 

-j;V»ndn Q-f"{u)udu 

"^^^^''''^ = -i::nu)dn = !::-f"{u)du = ^-z^^^'^^) 

JU2 J V / / N 

(3) We bound a from above: 

, , C/'»ndu «2 C /"(..) da 

"<•"• = f:;r(u)du < Qr'(u)du = ■ 

A similar argument bounds a from below. 

(4) We show that a{ui,U2) is strictly increasing in the second argument. 
Let ui < U2 < U3, Ui € [ii^,Uy]. Then 

Qr{u)udu+Qr{u)udu 



a{ui,U3) 



i{ui,U2) /„7 f'iu) du + a(n2, ng) Q f"{u) du 



Due to property ([3]) we have that 0(^1,^2) < U2 < 0(^2,^3). Thus 

a{m,U2)f:^ f'iu) du + a{m,U2)f:^ fin) du 
a(ni, na) > /"^ ///(^) dn " 0(^x1,^2) . 

A similar argument shows the result for the first argument. 
(5) Since ui < a{ui,U2) < U2 for all ui 7^ U2, we have (by the Sandwich 
Theorem) that 

u = lim ui < lim a{ui,U2) < lim U2 = u . 
Therefore, lim a{ui,U2)=u. □ 

Ul,U2 — 
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3. Interpolation and Particle Management 

The time evolution of equation ([1]) is described by the characteristic move- 
ment of the particles ([6j). Particle management is an "instantaneous" oper- 
ation (i.e. happening at constant time) that allows the method to continue 
stepping forward in time. It is designed to conserve area: The function value 
of an inserted or merged particle is chosen such that area is unchanged by the 
operation. A simple condition guarantees that the entropy does not increase. 
In addition, we define an interpolating function between two neighboring 
particles, so that the change of area under the interpolating curve satisfies 
relation (jlip . This interpolation is shown to be an analytical solution of the 
conservation law. 

3.1. Conservative Particle Management. Consider four neighboring par- 
ticles located at xi < X2 < X3 < xJ3 with associated function values ui, U2, 
U3, U4. Assume that the flux / is strictly convex or concave on the range 
of function values [minj(uj), maxj(?Xj)]. If U2 7^ M3, the particles' veloci- 
ties must differ /'(U2) 7^ fiu^), which gives rise to two possible cases that 
require particle management: 

• Inserting: The two particles deviate, i.e. f'{u2) < /'{us). If X3 — 
X2 > f^max for some predefined maximum distance dmaxj we insert 
a new particle {x23,U23) with X2 < X23 < X3, such that the area is 
preserved: 

{x23 - X2) a{u2,U23) + (a;3 - 3:23) a(^i23, ^^3) = (2:3 - X2) a{u2,U3) . (15) 
One can, for example, set X23 = ^2±M and find n23 by (fTSj) . or set 
U23 = and find X23 by 

• Merging: The two particles collide, i.e. /'(U2) > fiu's)- 1^x3— X2 < 
dmin for some predefined minimum distance (dmin = is possible), 
we replace both with a new particle (x23,U23) with X2 < 2:23 < X3, 
such that the area is preserved: 

(x23 - xi) a(ni, M23) + (xi - X23) a{u23,U4,) (16) 
= (x2-2;i)a(ni,M2) + (2:3 -2:2) 0(^2, ^3) + (X4-X3) 0(^3, U4) . 

We choose X23 = ^^i^, and then find M23 such that is satisfied. 

Figure [1] illustrates the merging step. 
Observe that inserting and merging are similar in nature. Conditions (|15|) 
and (|16p for ^23 are nonlinear (unless / is quadratic, see Rem. [20]) . For 
most cases U23 = "^"^"^ is a good initial guess, and the correct value can be 
obtained (up to the desired precision) by a few Newton iteration steps (or 
bisection, if the Newton iteration fails to converge). The next few claims 
attest that there is a unique value n23 that satisfies (fT5]) and (fT6l) . respec- 
tively. 



If more than two particles are at one position (x), all but the one with the smallest 
value (u) and the one with the largest value (u) are removed immediately. 
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Lemma 3. The function value U23 for the particle at X23 for equations (fTSll 
and (fTBl) is unique. 

Proof. From Lemma [2] we have that both a{ui,-) and a{-,U4) are strictly 
increasing. Thus, the LHS of both (|15|) and p6|) are strictly increasing (in 
U23), and cannot attain the same value for different values of U23. □ 

Lemma 4. There exists U23 € [u2,U3] which satisfies (jl5p . 

Proof. WLOG we assume U2 < U3. We define 

A = {X3 - X2)a{u2,U3) 

B{u) = {X23 - X2)a{u2,u) + (X3 - X23)a{u, U3) . 

So equation (fT5|) can be recast as B{u23) = A. The monotonicity of a(-, •) 
implies that B{u2) < A < B{u3). Since a is continuous, so is B, and the 
result follows from the Intermediate Value Theorem. The proof for U2 > U3 
is trivially similar. □ 

Lemma 5. If X2 = X3 = X23, there exists U23 G [ii2,'fi3] which satisfies ()16p . 

Proof. The proof is identical to the proof of Lemma [4] with the following 
definition of A and B{u): 

A = (x2 - xi) a(ui, M2) + {xi - X2) a{u3,Ui) 
B{u) = {x2 — xi) a(ui, u ) + (x4 — X2) a{u , U4) . 

□ 

Corollary 6. // particles are merged and inserted according to equations 
I j5i \10\i . then the total variation of the solution is either the same as before 
the operation, or smaller. 

Merging points only when X2 = X3 can be overly restrictive. The following 
theorem grants a little more freedom. 

Theorem 7. Consider four consecutive particles {xi,Ui) i = 1, . . . ,4. If 

1^*3 — ^^21 . /max I f"| \ ^ maxj Uj — mim Uj , „, 

^— > 4 . ; ' ' ' 17 

X3 — X2 \mm|/"|/ mm (a;4 — X3, X2 — xij 

then merging particles 2 and 3 with X23 = 22±£3 yi^ldg y^^^ ^ [^2, U3]. 

The min and max of /" are taken over the maximum range of ui, . . . , M4. 
Condition (jl7l) is trivially satisfied if X2 = X3. 

The idea of the proof is to consider merging in two steps. First, we find a 
value u such that setting U2 = U3 = u (while leaving X2 and 2:3 unchanged) 
preserves the area. Next, we merge the two particles to one with value ^23 
located at X23. To prove the theorem we use two lemmas: Lemma [8] bounds 
u away from U2 and U3 (but inside [1*2, 1^3]) • Lemma [9] bounds |ii23 — u\ from 
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above. We define three "area functions" : 

A = (X2 - Xi)a{ui,U2) + (X3 - X2)a{u2,U3) + (X4 - X3)a{u3,Ui) , 

B{u) = {x2 - xi)a{ui,u) + (x3 - X2)a{u, u) + {x^ - X3)a{u, U4) , 

C(n) = {x2 - xi)a{ui,u) + (X3 - X2)| [a{ui,u) + a{u, U4)] + (X4 - X3)a{u, U4) . 

Here A is the area before the merge that needs to be preserved, B(u) is the 
area when the particles 2,3 have the value u, and C{u) is the area when 
particles 2, 3 have been merged to a single particle at ^^"^^^ with value u. 

Lemma 8. The value u for which B[u) = A satisfies u E [^2,^3] and 

min/"(n) ^ 



. 1 min(x3 - xi,X4 - X2) b3 - ""21 

Ui\ > 7; 

2 X4 — xi 



max/"(ti) 

Lemma 9. The value U23 for which C{u23) = A satisfies 
[x3 - X2) [max(ni, U2, U3, U4) - min(n2,U3). 

Xi — Xi 



fori = 2,3. 



\U23 - "ul < 2 



max f"{u) 



min f"{u) 



2 



The proofs of the last two lemmas are tedious and uninspiring; they are 
relegated to the appendix for the interested reader's perusal. 

of Thm. Starting from the hypothesis of the theorem, we find 

1 min (x4 — X3, X2 — Xi) \u3 — U2I / min \ f"\ ^ ^ 

2 Xi — xi \max \ f"\ 
max |/"| ^ ^ (maxj Ui - min^ Ui){x3 - X2) 



min \f"\ ) X4 — xi 

> 2 maxj/^i y (maxj-uj - min {u2,U3)){x3 - X2) 
~ \min|/"|y X4 — xi 

Using Lemmas [8] and [9l we obtain 

\u — Ui\ > \U23 — u\ , 

for i = 2,3. Since we also have that u € [u2,U3] (from Lemma [9|), we 
conclude that U23 S [n2,n3]. □ 

Remark 10. Due to Thm. \7\ the merging step is robust with respect to 
small deviations in the distance of the merged particles. This also holds for 
the case X2 > X3, given the distance 1x3 — X2I is sufficiently small. 

Theorem 11. The particle method can run to arbitrary times. 

Proof. Let = miujUj, u,j = maxjUj, and C = maxf^^^^^] \ f"{u)\ • (u^ — 
u^). For any two particles, one has — f'{ui)\ < C. Define Axi = 

Xj+i — Xi. After each particle management, the next time increment (as 
defined in Sect. II. ip is at least Ats > SlSi^Si^ jf -^g not insert particles, 
then in each merge one particle is removed. Hence, a time evolution beyond 
any given time is possible, since the increments Ats will increase eventually. 
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Before 

After 




'(3:3,^3) 



Figure 1. Merging two parti- 
cles yields a new particle with 
function value chosen "conserva- 
tively" . This implies that the area 
of the two "triangles" is the same. 



Afi/'(«,) 




AtJ'iu) 



Xi x(u) X-i Xsh 



Figure 2. To find the x— value of 
a particle with given n— value, one 
locates the shock and then travels 
to the current time with velocity 
f{u) as given by (|20l) . 



When a particle is inserted (whenever two points are at a distance more than 
f^max)! the created distances are at least ^jy^, preserving a lower bound on 
the following time increment. □ 

Remark 12 (Choice of distance parameters). If possible, the minimum 
particle distance should be chosen dmin = 0. However, a small positive value 
of dmin also leads to a working method. This generalization is of interest if 
the characteristic equations have to be solved numerically, such as in Sect. [71 
The maximum particle distance dmax gives the minimal local resolution 
of the method. If the initial data is sampled with a resolution of /i, then a 
good choice for the maximum distance is /i < dmax < 2/i. Here, the upper 
bound comes from the fact that after an insertion the local particle distance 
is halved. Throughout our numerical simulations we use dmax — \^ since 
this gives, on average, a distance of h between particles in a rarefaction. 
This is an estimate that results from solving \ (dmax + |dmax) = ^ for dmax- 

3.2. Conservative Interpolation. Expression (113p defines an area be- 
tween any two points. We show that this defines an interpolating function 
v{x) between the two points. While an interpolation is not required for 
the particle management itself, it is useful for plotting the numerical solu- 
tion, interpreting its properties, and including source terms. As derived in 
Sect. 12. It the area between two points (xi,ui) and (x2,ti2) equals 
rx-i 

/ f (x) dx = (x2 - a;i)a(ni, U2) , 

J X\ 

assuming that / is strictly convex or concave in [7^1,^2]. We define the 
interpolation by the principle that any point (x, v) on the function v = v{x) 
must yield the same area when the interval is split: 

(x - xi)a(tii, w) + (x2 - x)a{y, U2) = (x2 - xi)a{ui,U2) ■ (18) 
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If Ml = U2, the interpolant is a constant function. Otherwise, (fTHI) can be 
rearranged to yield 

X -xi ^ a{ui,U2) - a{u2,v) ^ f'jv) - f'jui) ^^^^ 
X2-X1 a{ui,v) - a{u2,v) /'(■U2) - /'(ui) ' 
The last equality in (jl9p follows from 

Lemma 13. For any ui,U2, U3 with ui 7^ U2, the nonlinear average satisfies 

aiui,U2) - a{u2,U3) ^ /(^s) - f'jui) 
a{ui,U3) - a{u2,us) f'{u2) - f'{ui) ' 

Proof. By definition of the average function, we have 

= a(ni,M2) [/'(«)] "J +a(n2,U3) [/'(u)] + a(M3, ui) [/'(«)] 

"2 I / \ / r \1 "3 r j^'r Al^zA / \ r /.// xi^is 



= a(ni,n2) [/'(^)] + a(n2, ns) [[f{u)]ll - [f'{u)]ll)-a{m,us) [f'{u)]ll 
= {a{ui,U2) - a{u2,U3)) - (a(ni,U3) - 0(^2,^3)) [f'{u)]'^^^ . 

Rearranging the terms proves the claim. □ 

Observe that condition (jlSp is identical to the condition for particle in- 
sertion (jlSp . which means that any newly inserted particle must be placed 
on the interpolant. Relation (jl9p defines the interpolant as a function x{v). 
This is in fact an advantage, since at a discontinuity xi = X2, characteristics 
for all intermediate values v are defined. Thus, rarefaction fans arise natu- 
rally if f'{ui) < f'{u2)- If / has no inflection points between ui and U2, the 
inverse function v{x) exists. For plotting purposes we plot x{v) instead of 
inverting the function. 

The interpolation (jl9p can also be derived as a similarity solution of the 
conservation law (l2|), as follows. If ui = U2, we define v{x) = ui. Otherwise, 
one has f'{ui) 7^ f'{u2)- As derived in Sect.[2l the solution either came from 
a discontinuity (i.e. it is a rarefaction wave) or it will become a shock (i.e. it 
is a compression wave). The time Ati until this discontinuity happens is 
given by ([3]). At time t + Ati the particles have the same position xi = X2 = 
Xsh, as shown in Fig. [2j At this time the interpolation must be a straight 
line connecting the two particles, representing a discontinuity at Xsh- We 
require any particle of the interpolating function {x,v{x)) to move with its 
characteristic velocity f'{v{x)) in the time between t and t + Ati. This 
defines the interpolation uniquely as 

xiv) = XI - At^ if'iv) - /'(m)) = XI + l^^^J^S^^{x2 - xi) , (20) 



which equals expression (|19p . 

Not only is this interpolation compatible with the evolution of area under 
a conservation law, it also is a solution: 

Lemma 14. Together with the characteristic motion of the nodes, interpo- 
lation ()20p is a solution of the conservation law ([5]). 
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Proof. Using that Xi{t) = f'{ui) for i = 1,2 one obtains 




Thus every point on the interpolation v{x, t) satisfies the characteristic equa- 



CoroUary 15 (exact solution property). Consider characteristic particles 
with xi{t) < X2{t) < ■ ■ ■ < Xn{t). At any time consider the function defined 
by the interpolation (j20p . This function is a classical (i.e. continuous) solu- 
tion to the conservation law ([5]). In particular, it satisfies the conservation 
properties given in Sect. 

This corollary breaks down when shocks occur. 

Theorem 16 (TVD). With the assumptions of Thm.^ the particle method 
is total variation diminishing, thus it does not create spurious oscillations. 

Proof. Due to Cor. 1151 the characteristic movement yields an exact classical 
solution, thus the total variation is constant. Particle insertion simply refines 
the interpolation, thus preserves the total variation. Due to Thm.[71 merging 
yields a new particle with a function value U23 between the values of the 
removed particles. Thus, the total variation is the same as before the merge 
or smaller. □ 

Remark 17. The particle method approximates the solution locally by 
similarity solutions, very similar to front tracking by Holden et al. [12]. Front 
tracking uses shocks (after approximating the flux function by a piecewise 
linear, and the solution by a piecewise constant function). In comparison, 
our method uses wave solutions, i.e. rarefactions and compressions. 

3.3. Computational Aspects. 

Remark 18 (Shock location). The particle method does not track shocks. 
Still, shocks can be located. Whenever particles are merged, the new particle 
can be marked as a shock particle. Thus, any shock stretches over three 
particles {xi,ui), {x2,U2), (x^jUs), with the shock particle in the middle. 
Before plotting or interpreting the solution, a postprocessing step can be 
performed: The shock particle is replaced by a discontinuity, represented 
by two particles (x2, ui), (ic2, U3), with their position X2 chosen, such that 
area is conserved exactly. This step is harmless, since an immediate particle 
merge would recover the original configuration. As Fig. [3] illustrates, this 
postprocessing locates the shock with second order accuracy. The numerical 
results in Sect. 18.11 second this. 



tion ([6]). 



□ 



Remark 19 (Shock Speed). Since the particle method locates shocks with 
second order accuracy (Lemma fTSl) , the average shock speed is recovered 
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Figure 3. The location of the reconstructed shock 
is second order accurate in h. The square particle is 
a "shock particle" 

with the same accuracy. However, the instantaneous speed of a recon- 
structed shock is only a first order accurate approximation to the true shock 
speed. The reason is that the shock position is chosen according to a locally 
constant solution. This yields an 0{h) error in the local shock height, and 
thus an 0{h) error in the shock speed, according to the Rankine-Hugoniot 
condition [7\. In addition, at each particle management, the position of the 
reconstructed shock jumps a distance of order 0(/i^). Note that Riemann 
problems are solved exactly by our method. 

Remark 20 (Quadratic flux function). The method is particularly efficient 
for quadratic flux functions. In this case the interpolation (j20p between 
two points is a straight line, and the average (jl4p is the arithmetic mean 
a{ui,U2) = "^"^"^ . Thus, the function values for particle management can 
be computed explicitly. 

Remark 21 (Computational Cost). An interesting aspect arises in the com- 
putational cost of the method, when counting evaluations of the flux func- 
tion / and its derivatives /', /". Particle movement does not involve any 
evaluations, since the characteristic velocity of each particle f'{ui) does not 
change. Consider a solution on t G [0, 1] with a bounded number of shocks, 
to be approximated by 0{n) particles. Every particle merge and insertion 
requires 0(1) evaluations. After each time increment (jH), 0(1) management 
steps are required. The total number of time increment steps is 0(n). Thus, 
the total cost is 0{n) evaluations, as opposed to O(n^) evaluations for Go- 
dunov methods. Note that this aspect is only apparent if evaluations of /', 
/" are expensive, since the total number of operations is still O(n^). 

4. Sampling of the Initial Data 

When no source terms are present, the method has two sources of error: 
the initial sampling, and the merging of particles which contributes to the 
error in the neighborhood of shocks. Since, after the initial sampling and 
away from shocks, the solution is evolved exactly, it is natural to look for 
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ways to reduce the error due to the initial sampling. In some applications, 
the initial function uq may be representable exactly by the interpolation 
poll . In other cases, it has to be approximated. In Sect. 14.11 we show 
how well no can be approximated, as more and more particles are used. In 
Sect. 14. 2^ we outline a strategy of initializing a given number of particles in 
order to obtain a good approximation. 

4.1. Error Convergence. 

Lemma 22. Consider a smooth function w{x) on an interval x G 
with \f"{iD)\ >C> Vui G — |^), . Let v{x) denote the interpolant 

()20|) between and ^. Then \v{x) — w{x)\ = 0{h'^), and X^02 l^(^) ~ 
w{x)\dx = 0{h^) i.e. the approximation is second-order accurate in both 
L°° and norm^. 

Proof. Substituting a Taylor expansion w{x) = wq + w'qX + ^Wqx"^ + 0{x^) 
into /' yields 

f'{w{x)) = f'iwo) + f"{wo)w'oX + i {f"{wo)w'^ + f"'{wo)w'o') x' + 0{x') . 

(21) 

Using (j2ip in (j20p yields for the interpolation v{x) 

f'{v{x)) =/'(u;(-|)) + (/'(u;(|)) - /'(u;(-|))) ^ 

=/'(^o) - \f"{w,)w',h + \ {f"{wM + f"\w^X) + Oik') 
+ {f"(wo)w',h + Oih')){f^ + l) 

=nwo) + nwo)w'oX + I {f"{wo)w'^ + f"'iwo)w'i) + Oih^) . 

(22) 

Comparing (f2T]l and (f22|) yields 

nw{x)) - f'ivix)) = \ {f"{wM + n^oWi) (x' - Ih'^ + o{h') . 

(23) 

Using the Mean Value Theorem, we obtain the estimate 

w{x)-v{x) = - — — '- + 0{h-') 24 

2 f"{w) 

for some function value iv between v and w (or w{—^) and w{^)). Since the 
denominator is bounded from below by C, we have our result for the L°° 
norm. The error of the interpolation in the interval satisfies 

f-^o P = 12 n^) • 

Thus completing the proof. □ 



The power 3 in the order of the integral is needed so that the global L error, which 
results from adding up the errors from all the intervals, is second-order. 
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The formula for the error "density" is used for our adaptive sampling in 
the next section. We therefore name this density e: 



Non-smooth functions can be approximated as well, if the discontinuities 
are known: 

Theorem 23. Consider a piecewise smooth function no(x) with finitely 
many discontinuities (at known locations). Assume further that \ f"{uo{x))\ > 
C > Vx. Then uq can be approximated with second-order accuracy, using 
the interpolations (pOj) . 

Proof. First, represent each discontinuity in uq exactly, using two particles. 
This consumes only a finite number of points, thus the asymptotic behavior 
is not influenced. Second, place the remaining particles equidistantly at 
{xi,uo{xi)). Since the jumps are represented exactly, the maximum error is 
second order by Lemma [22l □ 

For non-convex flux functions, presented in Sect. [U the flux function can 
have inflection points at particles. In an interval bounded by an inflection 
particle, the second order accuracy is, in general, lost. First order accuracy 
is guaranteed by the following 

Lemma 24. Consider a smooth function w{x) on an interval x € [0,h], 
with \ f"{id)\ > G {w{0),w{h)) (so that the interpolation (I20p exists). 
Then the interpolation v{x) between and h, given by ([20|) . is at least first 
order accurate. 

Proof. The interpolation v{x) is monotonous, hence one can bound 



Remark 25. If the initial function is such that the flux crosses an inflection 
point, the error attains the form Hf — 1^11^1 ~ aK^ — b\og{h)h^ as h ^ 0. The 
local error in a single interval abutting the inflection point is of order 



Hence, the approximation is less than second order, but greater than first 
order. 

4.2. Adaptive Sampling. Due to Thm.[23t the initial data can be approx- 
imated with second order accuracy using equidistantly spaced points. Yet, 
for a fixed number of points, a non-equidistant spacing can yield a better ap- 
proximation. The presented particle method is designed for non-equidistant 
points. Hence, adaptive sampling strategies can be easily used. 




(26) 





error of order 0{h' 



2log(/l)). 
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To minimize the error for a given number of points we use a particle 
density proportional to e~2(x). For further information on this see [5j, for 
example. We define the integral 

E{x)= r e^{i)di. 
Jo 

and sample n + 1 points at positions 

Xi = E-' (^W^) • 

In the example presented in Fig. [6l this type of adaptive sampling is shown 
to reduce the initial error by a factor of about 2. 



5. Entropy 

We have argued in Sect. 13.21 that due to the constructed interpolation the 
particle method naturally distinguishes shocks from rarefaction fans. In this 
section, we show that the method in fact satisfies the entropy condition 

T]{u)t + q{u) < (27) 

for a convex entropy function rj, if the shocks are resolved well enough 
during a merge step. The following lemma considers the Kruzkov entropy 
pair r/fc(n) = \u — k\, qk{u) = sign(n — k){f{u) — f{k)). Holden et al. [HI 
Chap. 2.1] show that if ([27|) is satisfied for rjk, Qk (for all k), then it is satis- 
fied for any convex entropy function. Using the interpolation ()20p we show 
that the numerical solution obtained by the particle method satisfies this 
condition. 

Lemma 26 (entropy for merging). Let xi < X2 = x^ < X4 be the locations of 
four particles, with particles 2 and 3 to be merged, and f" > oH i.e. U2 > U3. 
If the value U23 resulting from the merge satisfies ui > U23 > u^, then the 
Kruzkov entropy J \u — k\ dx does not increase due to the merge. 

Proof. Consider the segment [xi,X4]. Let u{x) and u{x) denote the interpo- 
lating function before and after the merge, respectively. The interpolating 
function u is monotone in the value of its endpoints, thus u{x) < ii{x) for 
X G [x2,X4], and u{x) > u{x) for x E [xi, j;2]- The function 

I+{x) = 

can be used to write |x| = x + 2/+(— x). We identify two possible cases: 
k < U23 and k > U23. In the first case, k < U23, we write 

pX4 pX4 pX4 

/ |n— A;| dx = / {u — k)dx + 2 / I^{k — u)dx. 

J Xl J Xl J Xl 

■^For the case /" < 0, all following inequality signs must be reversed. 
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Due to the definition of u we have 

fXi rX2 rX4, 



rXA l'X2 /•X4 

/ (u-k) dx + 2 I+{k-u) dx + 2 L 

J Xl J X\ J X2 



\{k — u) dx. 



Since k < u on [xi,X2] and that is non-decreasing, we get 

> (u-k) dx + + 2 I+{k-u)dx. 

J Xl J X2 

Using A; < n on [xi, X2] we replace the zero with a diff'erent integral 

rx4 i'X2 i'X4 



rX4 pX2 rX4 

/ (u-k) dx + 2 I+{k-u)dx + 2 j I+{k-u)dx 

J Xl J XI J X2 

pX4 pX4 pX4 

/ {u — k) dx + 2 / — -u) dx = / \iL — k\dx. 

J Xl J XI J XI 



The other option we have is k > U23- In this case the proof is quite similar, 
but we start with |A: — 7i| instead, and remember that on [xi,X2] we have 
that u > u: 

X4 pX4 pX4 

\k — u\dx= / (k — u) dx + 2 / I^{u — k)dx 

Xl J Xl J Xl 

j-X4 I-X2 rX4 

= (k-u) dx + 2 I+{u-k) dx + 2 I+{u-k)dx 

J Xl J Xl J X2 

f-X4 /•X2 

> (k-u) dx + 2 I+{u-k) dx + 

J Xl J Xl 

j-X4 /•X2 rX4 

= (k-u) dx + 2 I+{u-k) dx + 2 I+{u-k)dx 

J Xl J Xl J X2 

pX4 pX4 pX4 

= / (k — u) dx + 2 / I^(u — k)dx= / \k — u\dx. 

J Xl J Xl J XI 

This ends the proof. □ 

The assumption of Lemma[26l implies that shocks must be reasonably well 
resolved before the points defining it are merged. It is satisfied if the points 
to the left and right of a shock points are not too far. The condition can be 
ensured by an entropy fix: A merge is rejected a posteriori if the resolution 
condition is not satisfied. Then, points are inserted near the shock, and the 
merge is re-attempted. 

Remark 27. With the entropy fix, a merge does not necessarily reduce 
the number of points. Based on numerical evidence, we conjecture that the 
statement of Thm. [TT]remains valid, although its proof cannot be transferred 
in a straightforward fashion. 

Theorem 28. The presented particle method yields entropy solutions. 



EXACTLY CONSERVATIVE PARTICLE METHOD FOR CONSERVATION LAWS 18 





1 


(1) 




















(3) 



Figure 4. Particle management around an inflection par- 
ticle (/"(ua) = 0) results in one of three possible configu- 
rations. Each configuration allows for more area under the 
function than the previous one. Here we see three archetypal 
particle configurations that result. 



Proof. During the characteristic movement of the points, the entropy is con- 
stant. This is due to Cor. 1151 which tells that the interpolation is a classical 
solution to the conservation law. Particle insertion does not change the in- 
terpolation, thus it does not change the entropy. Merging does not increase 
the entropy if the conditions of Lemma [261 are satisfied. □ 

6. NoN-CoNVEx Flux Functions 

So far we have only considered flux functions without inflection points 
(i.e. /" always has the same sign) on the range of function values. In this 
section, we generalize our method for flux functions / where /" has a fi- 
nite number of zero crossings n| < • • • < . Between two successive points 
u G [it*,M*_|_i] the flux function is either convex or concave. We impose the 
following requirement for any set of particles: Between any two particles 
for which /" has opposite sign, there must be an inflection particle {x,u*). 
Thus, between two neighboring particles, / never has an inflection point, and 
the fundamental ideas from the previous sections transfer. In particular, the 
characteristic movement of particles is unaffected, and the interpolation be- 
tween two particles remains uniquely defined by ()20p . It has infinite slope at 
inflection points (since f"{u*) = 0), but this is mostly harmless. However, 
two complications arise. First, every proof that relies on having a lower 
bound on /" does not transfer easily. Second, merging particles when an 
inflection particle is involved requires a special treatment. The standard ap- 
proach, as presented in Sect. 13. H removes two colliding points and replaces 
them with a point of a different function value. If an inflection particle is 
involved in a collision, points must be merged in a different way so that an 
inflection particle remains. 

We present one such special merge for dealing with a single inflection 
point (we do not consider here the interaction of two inflection points). Also, 
for simplicity, we consider a collision with identical point positions. Since 
the inflection particle must remain (although its position may change), we 
consider five neighboring particles and not four as before. Let {xi,Ui),i = 
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1, . . . , 5 be these particles so that X2 = 2:3, f"{u^) = 0, and (WLOG) > 
0, i.e. the inflection particle is the slowest. The other cases are simple 
symmetries of this situation. The special merge consists of three different 
attempts to find the new particle configuration. The first two attempts may 
fail to provide a solution, in which case the next is attempted. 

(1) Remove particle 2 and increase such that area is preserved. Ac- 
cept, if 2:3 is not increased beyond X4. 

(2) Remove particle 2, set = X4 and increase both such that area is 
preserved. Accept, if X3 and X4 are not increased beyond X5. 

(3) Remove particle 4, set X3 = X5 and lower U2 such that area is pre- 



Theorem 29. One of the three attempts listed above will succeed. 

Proof. Following from continuity and monotonicity of the average function 
a(-, •), the three steps provide a continuous, monotonous increase in area. In 
the first attempt, the smallest area is achieved with X3 unchanged. This area 
is necessarily smaller than the original area (since one can also get here by 
lowering U2 to M3). The area increases as X3 is increased. The configuration 
with x^ = X4 has the maximum area for the first attempt, and the minimum 
area for the second. Again, the area increases as X3 = X4 increase. The 
configuration with x^ = X4 = X5 has the maximum area for the second 
attempt, and the minimum for the third. Area increases as the new value 
of U2 increases, and achieves its maximum value for an unchanged U2. This 
area is necessarily larger than the original area. Consequently, one of the 
attempts must succeed. □ 

Remark 30. The resulting configuration may involve a new discontinuity 
(since X3 = X4 or X3 = X5). However, this is not a shock, but a rarefaction, 
since the particles will move away from each other. Consequently, these 
particles should not be merged. 

The five-point particle management guarantees that in each merging step 
one particle is removed, as used in Thm. [TTl In Sect. 18. 2^ numerical results 
on the Buckley-Leverett equation are presented. Since each of the three 
attempts covers a non-overlapping range of areas, the resulting configuration 
is independent of the order in which they are attempted. 



An important extension of the conservation law ([T]) is to allow a source 
term in the right hand side. This can be a function of x, t, the function 
value u, and in principle also of derivatives Ux, Uxx, etc. In the current work 
we consider the simple balance law 



served. 



7. Sources 




(28) 
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The method of characteristics ^ yields an evolution for each particle 



With sources, the equation ceases to have exact conservation properties. 
Consequently, the interpolation derived in Sect. [3^2] is no longer a solution. 
While in special cases more complicated interpolation functions could be 
defined (depending on both / and g), here we construct an approximate 
method that is more general. Assume that the advection dominates over 
the source, which is the case in many applications. Thus, the interpolation 
and particle management are based solely on the flux function /. 

The source g results in a vertical movement of the particles during their 
Lagrangian evolution. While in the absence of sources the next time of a 
particle merge can be computed a priori, now we solve the particle evolution 
()29|) numerically, for instance by an explicit Runge-Kutta scheme. Merging 
takes place when two particles are too close (see Rem. [TO]l . In Sect. 18.31 we 
present numerical results. 

Remark 31. The balance law (|28p is solved correctly at characteristic 
points. Particle management, however, is based on an "incorrect" interpo- 
lation, since the source is neglected for the definition of area. The numerical 
results in Sect. 18.31 indicate that this does not cause problems for merging 
particles. However, inserting particles into large gaps may lead to significant 
misplacements, when the source is "active". Thus, with sources, insertion 
should either be avoided completely, or particles be adaptively refined. We 
shall address the important aspect of adaptivity in future work. 

The presented approach incorporates sources directly into the character- 
istic equations. An alternative approach is operator splitting: First move 
particles neglecting the source, then correct function values according to 
the source. While the characteristic method is more precise, the splitting 
approach is more general. In particular, it can deal with source terms that 
involve derivatives of u. 



The presented particle method is applied to various examples. In all 
cases, the "exact" reference solution is obtained or verified by a high resolu- 
tion CLAWPACK p] computation. We compare the accuracy of the particle 
method with numerical solutions obtained by CLAWPACK, considering sim- 
ilar resolutions. By construction, the particle method does not keep a fixed 
resolution. To compare resolutions we use the same number of particles to 
initialize the particle method as the number of cells in the corresponding 
CLAWPACK run. By keeping dmax = we find that the number of parti- 
cles remains more-or-less constant throughout the computation. Shocks are 
located via post-processing before the error is measured. 




(29) 



8. Numerical Results 
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Figure 5. The particle method for /(u) = before and 
after shocks arise. The gray sohd hne is the hi-resolution 
solution from CLAWPACK and the dashed line is the initial 
condition. 



In Sect. 18. H the evolution and the formation of shocks of smooth initial 
data under a convex flux function are considered. The convergence error 
before and after the occurrence of shocks is investigated numerically. In 
Sect. 18.21 as an example of a non-convex flux function, the Buckley-Leverett 
equation is considered, and in Sect. 18. 3^ Burgers' equation with a source is 
simulated. The source code and all presented examples can be found on the 
particleclaw web page |19j . 

8.1. Convergence Error. Figure[5]shows the smooth initial function no(a;) = 
0.5 + 0.2exp(— x^) cos(7r2;), and its time evolution under the flux function 
f{u) = \u^. Initially, we sample points on the function uq. At time t = 1, 
the solution is still smooth, thus the particles lie exactly on the solution. 
By the time t = 10, a shock has emerged and interacted with a rarefaction. 
Although the numerical solution uses only a few points, it represents the 
true solution well. 

From this example, the numerical accuracy of the particle method is ex- 
tracted. For a sequence of particle densities, the initial data are sampled 
twice: equally spaced and adaptively. The particle method is applied with 
post-processing, as described in Rem. [181 The error is measured in the 
true norm for function, which is possible due to the interpolation ()20p . 
Figure [6] shows results for initial sampling error, and error after a time evo- 
lution. Initially (t = 0), the approximation is second order accurate for both 
sampling strategies (see Thm. [23]) . The advantage of the adaptive sampling 
is evident from the lower error that it creates in the interpolation. 

After shocks have occurred {t = 10), the approximate solution without 
locating shocks is only first order accurate, since at any shock an error of the 
order height x width of the shock is made. However, the post-processing step 
recovers the second order accuracy. Hence, the particle method is second 
order accurate, even at locating shocks. One also sees that the advantage 
gained initially from the adaptive sampling is nearly lost at t = 10. 
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Figure 6. Error of the particle method when solving for 
the flux function initially, and after a shock has been 
formed. The figures compare equidistant with adaptive sam- 
pling and show the error that would result from not sharp- 
ening the shock. 



For this example, CLAWPACK yields results of similar accuracy, as shown 
in Fig. [71 Since CLAWPACK is a method for calculating cell averages, we 
cannot find the true error. Instead, given a coarse-grid calculation and 
a fine-grid reference solution, we calculate the error in the area (function 
value times cell-size) for each of the coarse cells using the fine-grid solution. 
Adding all these errors together gives the relevant error. One can see 
that the CLAWPACK solution drops to first-order accuracy around shocks, 
which is due to numerical dissipation. To investigate the accuracy away 
from shocks, we also consider the error while ignoring a small fixed domain 
surrounding each shock. The same error measure is also applied to the error 
calculations of our particle method. Of course, since post-processing already 
yields second order accuracy, this only reduces the size of the error, and does 
not change the order of convergence, as it does with CLAWPACK. 

8.2. Non-Convex Flux Function. As an example of a non-convex flux 
function, we consider the Buckley-Leverett equation 

ut + {f{u))., = , with f{u) = u^/{u^ + i(l - uf) , (30) 

which is a simple model for two-phase fluid flow in a porous medium (see 
LeVeque jl5l). We consider piecewise constant initial data with a large 
downward jump crossing the inflection point, and a small upward jump. 
The large jump develops a shock at the bottom and a rarefaction at the top, 
the small jump is a pure rarefaction. Around t = 0.2, the two similarity 
solutions interact, thus lowering the separation point between shock and 
rarefaction. Figure [8] shows numerical results. The solution obtained by 
the particle method (dots) is compared to a second order CLAWPACK 
solution (circles) of similar resolution. The particle method captures the 
behavior of the solution better; in particular, the rarefaction is represented 
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Figure 7. A comparison of the errors given by the particle 
method and CLAWPACK when solving for the flux function 
ju"^. Without removing the errors from the shock region, 
CLAWPACK is only first-order accurate. 
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Figure 8. Numerical results for the Buckley-Leverett equa- 
tion at various times t 



very accurately. Only directly near the shock are inaccuracies visible. The 
solution away from the shock is nearly unaffected by the error at the shock. 

Numerical results show (see Fig. [9]) that both CLAWPACK and the pre- 
sented particle method do not achieve second order accuracy for this prob- 
lem. Nevertheless, the particle method has a much better accuracy than 
CLAWPACK. The drop in accuracy is, presumably, due to inflection point 
in the Buckley-Leverett flux function, similar to the drop in accuracy of the 
sampling outlined in Remark 1251 

8.3. Source Terms. We consider Burgers' equation with a source 

ut + {^u^), = b'{x)u. (31) 

It is a simple model for shallow water flow over a bottom profile b{x). As in 
|14j . we consider the domain x G [0, 10], and choose 




cos(7rx) X G [4.5, 5.5] 
otherwise . 
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Figure 9. A comparison of the errors given by the particle 
method and CLAWPACK when solving the Buckley-Leverett 
equation (|3Up . While the particle method gets significantly 
better results than CLAWPACK (with or without the errors 
from around the shock), both methods are less than second 
order accurate. 

The source term is included into the method of characteristics, as explained 
in Sect. [71 The time stepping is done by a fourth order Runge-Kutta scheme. 
Figure [TU] shows the computational results. The particle method (dots) ap- 
proximates the solution significantly better than the second order CLAW- 
PACK scheme (circles). A particular aspect in favor of the characteristic 
approach is the precise (up to the resolution of the ODE solver) recovery of 
the function values after the obstacle. Since particles are moved indepen- 
dently according to the characteristic equations, an accurate time integration 
obtains the function values after the obstacle almost exactly, independent 
of the resolution of particles. Note that an efficient treatment of the source 
requires a special consideration of its discontinuities, either in the quadra- 
ture of the source (finite volume) , or in the integration of the characteristic 
ODE (particle scheme). 



We have presented a particle method that combines the method of char- 
acteristic, local similarity solutions, and particle management to a numerical 
scheme for for ID scalar conservation laws. The method conserves area ex- 
actly. It is TVD, yet second order accurate, even at locating shocks. It 
performs promisingly in various examples, as the numerical comparisons 
with a second order finite volume scheme show. 

The particle method is an interesting alternative to fixed grid approaches, 
whenever conservation of mass is crucial, or shocks need to be located accu- 
rately. In addition, entropy is reduced only when particles are merged, which 
makes the approach suited for applications in which the evolution of mass 
and energy has to be reflected as precisely as possible. Furthermore, the 



9. Conclusions and Outlook 
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Figure 10. Time evolution of Burgers' equation with a 
source as given by ([3T|) . The dots, circles, gray line and 
dashed line are, respectively, our particle method, CLAW- 
PACK, the high-resolution solution and the initial conditions. 



method yields good results when few particles are used, in particular shocks 
between nearly-constant states are located well. This makes the approach 
attractive whenever scalar ID conservation laws arise as sub-problems in a 
large computation, and only a few degrees of freedom can be devoted to the 
numerical solution of a single sub-problem. Examples are flows in networks 
(e.g. car traffic), and PDE constrained optimization. 

As a first generalization, we have included source terms in the scheme. 
The method, still based on the method of characteristics, yields solutions 
of rather striking accuracy, compared to classical finite volume schemes. In 
future work, more general source terms will be considered, such as non-local 
convolutions, and terms involving derivatives of the solution. In these cases, 
the method of characteristics has to be replaced by a more general splitting 
approach. 
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Figure 11. A comparison of the errors given by the particle 
method and CLAWPACK when solving for Burgers' equation 
with a source term (|31|) . at t = 6 (after the two shocks have 
interacted). Here too, the particle method provides a better 
accuracy than CLAWPACK. 

Fundamental steps towards a more powerful particle method will be the 
generalization to higher space dimensions and to systems of conservation 
laws. Problems in multiple dimensions can be approximated by ID prob- 
lems using fractional steps. In this sense, the particle scheme could replace 
classical ID Riemann solvers by ID wave solvers. However, this approach 
is not fully satisfactory, since due to the required remeshing steps the ben- 
efits of a meshfree particle approach may be lost. On the other hand, with 
truly meshfree approaches in 2D/3D, one has to address the challenge that 
particles forming a shock need not necessarily collide. Possible remedies 
are the introduction of a numerical pressure, or the tracking of an unstruc- 
tured triangular mesh. The movement of particles according to the method 
of characteristics can also be interpreted as a moving mesh approach [Ij. 
Thus, ideas from this area could lead to particle strategies in higher space 
dimensions. 

With systems, one difficulty is the presence of multiple characteristic ve- 
locities. One approach is to choose one Lagrangian velocity, which need 
not be a characteristic velocity. Coupling terms that appear in the moving 
frame equations are treated as source terms for each individual equation. 
Alternative approaches may use exact similarity solutions of the full system 
as building blocks. In this case, a single set of particles may not suffice, 
since two neighboring similarity solutions may interact. 
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Appendix A. Appendix 
The proofs of Lemmas [8] and [9] use a short lemma: 

Lemma 32. The derivative of a{u,v) with respect to either of its variables 
is hounded from below and above as follows: 



1 / min/'' 



2 \ max /" 



< 



da \ da 



< 



1 / max/"^ ^ 

2 V mill /" 



Here max/" and min/" are taken over the interval \u, v\. 
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Proof. This lemma follows from the definition of a: 

da, ^_nu)f:f"iu;)iu;-u)du; 



du 



{u,v) 



2 



^ max /" \ oj — u du) ^ I I' max /" 



min/"y (w — u)^ 2 ymin/" 
The other bounds (on ^ and the lower bound) have similar proofs. □ 

of Lemma\^ WLOG we assume that > U2, and show for U2- We bound 
A — B{u2) from below: 

A - B{u2) = (x3 - X2)(a(u2, M3) - a{u2, U2)) + {x^ - X3){a{u3,U4,) - a{u2,Ui)) 

> (X3 - X2){us - U2) min ^(n, v) + {x4^ - X3)(n3 - U2) min ^(n, 

ov ou 

N 1 / niin f " \ ^ , , 

>(,._,,,(„3_„,)_(^_I_j . (32) 

Since we are looking for u such that B{u) = ^, the previous bound is also a 
bound on B{u) — B{u2)- From the Mean Value Theorem we have ^ € [u, U2] 
for which 

B{u) - B{u2) 
B'iO 

B{u) - B{u2) 



{X2 -Xl)^{i,Ui) + (X4 - X3)|^(C,n4) + (X3 " X2) 

^ B{u) - B{u2) 

"(^4-x0(^)' 

In the last step we used the upper bound on |^ and that '^^j,, > 1- From 
(f32|) we conclude that 

1 (X4 - X2)(ti3 - W2) / min/" ^ ^ 

U — U2 > -- 



2 X4 — xi \max /" 

Similarly, one can show that U3 - u > 1 (^3-xi)(u3-»2) / EliRfl)\ □ 

of Lemma Again, WLOG we assume that > U2- This time we first 
bound \C{u) — A\ from above: 

\C{u)-A\ = C{u)-B{u) 

= -^—^ — ^(a(ui, -u) + a({t, U4) — 2a(u, u)) 
< (x3 - X2) [max(ui) - min(n2, n3)] . 
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Recall that C(u23) = A, and that C > due to the monotonicity of a. 
Thus, for some point ^ 

|C(u)-C(t.23)| 



\U - U23\ 



^ (x3 - X2) [max(n^) - min(n2,^3)] 
~ min C 

^ ^ (x3 - X2) [max(ttj) - min(tt2,n3)] / max/^^ ^ 
~ {X4-X1) \mmf" 



□ 



